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Abstract — Extensive simulations of the ultra-high energy cos- 
mic ray (UHECR) propagation in the Galactic magnetic field 
(GMF) have been performed, and the results are presented. 
The use of different available models of the large-scale GMF 
and/or primary particle assumptions leads to distinctly different 
deflection patterns of the highest energy cosmic rays (CR). The 
lensing effects of the Galactic field modify the exposure of 
an UHECR experiment to the extragalactic sky. To quantify 
these effects for the Pierre Auger experiment, we performed a 
correlation analysis of the simulated cosmic ray event samples, 
backtracked from the Earth to the Galactic border, with the 
active galactic nuclei (AGN) from the 12th edition of the Veron- 
Cetty and Veron catalogue. Further forward-tracking studies 
under plausible UHECR sources scenarios are needed to allow 
for direct comparison with the observed correlation between the 
nearby AGN and the highest energy Auger events. 

I. Introduction 

Magnetic fields in the Milky Way and other galaxies are 
investigated by means of measurements of the Faraday rotation 
of the polarized light from pulsars and extragalactic sources, 
and through detection of synchrotron radiation emitted by rel- 
ativistic electrons [1], [2]. The observations of spiral galaxies 
like ours suggest that there is a large-scale field, that in the first 
approximation follows the spiral arm structure, and is parallel 
to the galactic disk. In the vicinity of the Sun such large-scale 
field is pointing approximately in direction of galactic latitude 
i = 80°. The spiral field extends above and below the galactic 
disks on a kiloparsec scale and forms a kind of halo, which is 
confirmed by observations of the galaxies that are seen edge- 
on. A random Galactic magnetic field, of the strength similar 
to that of the regular component, has also been observed [1]. 
The coherence length of the turbulent field is at most of the 
order of several tens of parsecs. The average strength of the 
total magnetic field is about 6 jiG near the Sun, and increases 
towards the Galactic Center region. The observations of the 
Galactic diffuse soft X-rays indicate on a possible CR and gas 
pressure driven wind, which would transport magnetic fields 
from the disk to the halo [3], [4]. 

The observed Faraday rotation measures (RM) allow to 
probe the Galactic magnetic field in the Solar proximity [5], 
[6], [7]. However, the limited data sample makes the 
interpretation of these observations difficult. In addition, the 
RM measurements provide the line-of- sight convolution of the 
magnetic field with the thermal electron density, which in turn 
is not well known. New radio polarization observations, with 
better sensitivity and angular resolution are needed to obtain 



more conclusive picture of the magnetic field distribution in 
the Milky Way. A significant improvement of the data will 
be provided by new-generation radio telescopes. One of such 
instruments, the Low Frequency Array (LOFAR), is already 
under construction [8]. LOFAR will be in particular able to 
trace the magnetic fields in halo regions of the Milky Way 
and other galaxies, through detection of a few meters radio 
synchrotron emission from low energy cosmic rays [9]. 

The Pierre Auger Observatory [10] provides a new and inde- 
pendent way of studying cosmic magnetic fields, by collecting 
cosmic ray events at the extreme energies (above lOEeV = 
10 19 eV) with unprecedented statistics and data quality. Re- 
cently, the Pierre Auger Collaboration observed [11], [12] 
a significant anisotropy in the arrival directions of cosmic 
rays above ~ 60EeV. These cosmic rays correlate over 
angular scales less than 6° with the locations of nearby (D < 
100 Mpc) AGN from the 12th edition of the Veron-Cetty and 
Veron (VCV) catalogue [13]. The Pierre Auger experiment 
has also detected [14] a strong steepening of the cosmic ray 
flux above 4 x 10 19 eV. Both observations are consistent with 
the standard scenarios of the UHECR production in the extra- 
galactic astrophysical acceleration sites and thus represent the 
important step towards the "charged particle astronomy". 

The final way the CR astronomy will be done will strongly 
depend on the primary mass composition of the highest 
energy cosmic rays, since the magnetic deflections scale in 
proportion to the CR atomic number Z. The elaboration 
of relevant analysis methods requires detailed investigation 
of the UHECR propagation in cosmic magnetic fields, and 
first of all in the Galactic field. The different aspects of 
the UHECR propagation in the GMF have been extensively 
covered in the literature, the non-exhaustive list of previous 
results can be found in [15], [16], [17], [18], [19], [20], 
[21]. We present in this paper results of our own studies, 
performed in the light of the observed AGN correlation. We 
used the standard method of CR backtracking (see Sec. [II]). 
Three distinctive large-scale field models have been chosen 
(Sec. [III]). A large number of CR events has been simulated 
using energy spectrum, arrival direction distribution, and mass 



composition described in Sec. IV and propagated under the 
assumed GMF models. Resulting magnetic deflections are 
presented in Sec.JV] Modification of the extragalactic exposure 
due to the large-scale GMF lensing effects is discussed in 



Sec. VI followed by conclusions. 
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Fig. 1. The spiral disk field in the HMR models [16]. The volume of galactic 
region with magnitude of magnetic field larger than 1 fiG is shown in yellow- 
green. The position of the Solar system is denoted by a red sphere. 



II. CR Backtracking in the Galactic field 

Since the energy loss length in the local Universe largely 
exceeds the size of the Galaxy (see, for example, [22], [23], 
[24]), we can neglect the energy losses when considering the 
propagation of the ultra-high energy cosmic rays through the 
Milky Way. This makes it possible to inverse the problem 
and, instead of following the trajectory of a particle from 
the Galactic halo border till the arrival on Earth from a 
certain direction, to propagate the corresponding anti-particle 
from Earth in the same direction (backtracking method). The 
trajectory of a UHE cosmic ray is then obtained by integration 
of the equations of motion of an ultra-relativistic particle with 
the energy E and charge q in a quasi- static magnetic field B\ 



dr 
dct 



dv 
dct 



qc 



v x B , 



(1) 



where v = p/\p\ , p = Edr/c 2 dt, is the unit vector parallel to 
the velocity of the cosmic ray. Since in the highly structured 
Galactic magnetic field the field strength can vary significantly 
along the cosmic ray trajectory, we have implemented for the 
integration of the equations of motion ([]]) the Runge-Kutta 5 th 
order scheme (RK5) with the adaptive step size control [25]. 
To avoid the "numerical dissipation" of the cosmic ray energy 
E (the effect studied in detail in [26]), we imposed the energy 
conservation during the propagation steps, by preserving the 
absolute value of the particle's velocity vector. 

The integration accuracy, represented by the accepted trun- 
cation error in the RK5 technique, has been optimized correlat- 
ing the backtracked directions with the selected astrophysical 
objects. An accuracy level of 10 ~ 6 has been adopted, for 
which the changes in the backtracked directions become 
negligible with respect to the chosen step in angular separation 
cUax from the catalogue objects (see the section fVI-A[). 



The cosmic ray trajectory has been followed till the galac- 
tocentric distance of 20 kpc (the Galaxy "border"), beyond 
which the field strength is supposed to be negligibly small. 



III. Considered large-scale GMF models 

A strong large-scale magnetic field with a stationary or 
oscillating configuration can be generated from a weak seed 
field as a result of the non-uniform (differential) rotation of 
the galactic gas with strong turbulent motions [2] . The regular 
field structure obtained by this so-called dynamo mechanism 
can be described by modes of different azimuthal symmetry 
in the disk, and vertical symmetry perpendicular to the disk 
plane. The strongest mode is the spiral disk field with either 
7r (bisymmetric, or BSS) or 2tt (axisymmetric, or ASS) sym- 
metry. Other weaker field modes can also be excited. As for 
the vertical symmetry, the large-scale galactic fields can either 
keep the direction while traversing the disk plane (be of even 
parity), either change it to the opposite (odd parity modes). 
Both symmetry properties are reflected in the toponym of the 
large-scale field models. Combined with the spiral structure 
type, 4 possible spiral field patterns are denoted BSS-S, BSS- 
A, ASS-S, ASS-A, where -S (symmetric) stands for the even 
parity, and -A (antisymmetric) - for the odd parity. 

We have chosen amongst many available large-scale GMF 
models three typical ones with distinct qualitative differences. 
Two models are the spiral disk field models of bisymmetric 
even parity, and axisymmetric odd parity from the paper [16] 
by Harari, Mollerach and Roulet (HMR). These models (see 
Fig.[T]) represent a version of the Stanev model [15], smoothed 
out in order to avoid the field discontinuities. 

The third model is a modification of the model [19] by 
Prouza and Smida (PS), that has been proposed by Kachel- 
rieB et al. [20]. In addition to the bisymmetric even parity 
spiral field, of the structure similar to the one in the HMR 
BSS_S model, it features two additional large-scale halo GMF 
components: toroidal azimuthal field above and below the 
Galactic disk, and poloidal (dipole) field. 

IV. Parameters adopted for simulated CR events 

We simulated a large number (stated in Tab. [I] and |TTJ> of 
cosmic ray events with the parameters below. 

A. Energy spectrum 

The energy values of the simulated events were bounded 
between 40EeV and 150EeV, and distributed according to the 
Auger measurement [14], as a power law E~ 4 - 2 . 

B. Angular distribution 

Two options for the angular distribution of the events have 
been considered: 1) Uniform distribution over the whole sky, 
and 2) Distribution according to the exposure of the Auger 
Surface Detector (SD) for zenith angles 6 Z < 60°. For the 
considered energies the Auger SD acceptance area is saturated, 
which produces a simple analytic dependence of the exposure 
upon declination [27]. 

C. Primary mass composition 

4 pure compositions have been considered: protons (Z = 1), 
carbon (Z = 6), silicon (Z = 14), and iron (Z = 26) nuclei. 



(a) Bisymmetric even parity spiral field 



(b) Axisymmetric odd parity spiral field 



Fig. 2. Cosmic ray deflection patterns expected for protons under HMR models [16]. The deflections shown represent the arcs of the great circle joining the 
simulated arrival direction on Earth, and the backtracked one at the Galaxy border (indicated by an arrow head). To visualize deflections, every 100^ event 
has been drawn from the large simulated event data sample described in the section [TV] The dashed line denotes the supergalactic plane. The thin solid lines 
show the equal integrated Auger SD exposure sky regions within the detector field of view for zenith angles Z < 60°, delimited by the black solid line. 




V. Magnetic deflections of cosmic rays 

The cosmic ray deflection patterns expected from the three 
selected GMF models are very different, as it can be seen from 
Figs. [2] and [3] on the example of primary protons. Except for 
the directions towards the Galactic Center, where the character 
of deflections is complex due to the multiple field reversals, 
there are trends, typical for each model. In the bisymmetric 
even parity spiral disk field, there is a flow in the direction of 
the Galactic North Pole, i.e. cosmic rays in general arrive at 
the Galactic halo border from higher Galactic latitudes than 
observed on Earth. In the odd parity field, the backtracked 
arrival directions in each Galactic hemisphere are shifted to 
the poles. This field configuration would obscure the Earth 
sky from the cosmic rays coming from the directions close to 
the Galactic plane. The presence of dipole and toroidal fields 
makes the deflection pattern even more complex (see Fig. [3]). 

The cumulative distributions (c.d.f.) of resulting deflections 
for the HMR BSS_S model are shown on Fig. [4] The deflection 
values for the three considered GMF models are summarized 
in the table [TJ by means of percentiles at 50% (median), 
68.27%, and 95.45% of the c.d.f. Neglecting the distribution 
tails, one can see that the deflection values scale rather well 
with the atomic number Z of primary nuclei. Despite the 



Fig. 4. C.d.f. of the deflection angle values (in °), expected in the HMR 
BSS_S model [16] for the four assumed mass compositions. 
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1 
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4.3 
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75.9 
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HMR axisymmetric odd parity model 
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100000 


5.4 


6.6 


10.9 
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100000 


31.5 


40.2 


77.9 


Si 
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100000 
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97.6 


139.2 


Fe 


26 


100000 


84.0 


100.4 
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PS model version by Kachelriefi et ah 



p 


1 


20000 


3.0 


4.3 


20.8 


c 


6 


20000 


16.7 


22.9 


64.7 


Si 


14 


20000 


37.6 


47.2 


85.7 


Fe 


26 


20000 


57.5 


70.9 


118.5 



TABLE I 



Magnetic deflections (in °) at the indicated percentiles of the 
c.d.f. (shown on flg|4]for one of the three gmf models). 



additional halo fields in the PS model version, the higher 
overall normalization of the disk field strength in the HMR 
models, and its more important halo extension [16], [19], [20] 
are the reasons why in average the deflections for the HMR 
models are larger. 

It is instructive to compare the obtained UHECR deflections 
(see Fig. [4} with those expected for turbulent field with a 
coherence length / CO h. In such a case, on a distance d ^> / CO h 
a cosmic ray will propagate in the random walk regime, with 
the resulting r.m.s. deflection angle given by [22], [23] 



proton 







random \ 



V^TkH V%oh [loo pc] B \pG\ 
' L13 Z £[10 EeV] ' (2) 

where the effect of dynamical friction [22] is taken into ac- 
count. The path d can be roughly estimated using the simulated 
events (see Sec.[iV|). For the HMR BSS_S field configuration, 
the median distance <i 50 % traveled till the Galaxy border 
is 21.0 kpc for protons, and 28.3 kpc for iron nuclei, while 
the respective values for the two primaries corresponding to 
95.45% of the c.d.f. are 28.1 kpc and 70.3 kpc. Under reason- 
able assumption of the equal field strength for the large-scale 
and turbulent components, an r.m.s. turbulent field strength 
is expected to decrease with a distance to the Galactic plane 
at a kiloparsec scale, which is much smaller than the above 
mentioned d values. Hence, for the energies of interest and for 
most of the directions in the sky the UHECR deflections from 
the random Galactic field component ^ will be considerably 
smaller than the ones from the large-scale GMF (see also [28]). 
Therefore, the former can be neglected at the first approach, 
though a realistic simulation of the UHECR propagation in the 
Galactic magnetic field has to take into account the turbulent 
field. 

VI. GMF EFFECTS ON THE EXTRAGALACTIC EXPOSURE 

In addition to magnetic deflections, (de-)focusing effects of 
the Galactic magnetic field will modify the exposure of an 
experiment to the extragalactic sky [16], [17], [20]. Though 
the Liouville theorem holds and the isotropic cosmic ray 
flux outside the Galactic halo will remain isotropic on Earth, 
the correspondence between the arrival directions on Earth, 
and the sources contributing to the flux on the halo border 
will be quite complex. To study how this correspondence is 
realized, we have used the backtracked simulated events from 



the section [IV] We mapped the decimal logarithm of the ratio 
of the number of events on the halo border (r = 20 kpc) to the 
one on Earth, using the HEALPix equal area celestial sphere 
pixelization [29]. These maps, shown on Figures [5j [6] and [7] 
allow to estimate the (de-)magnification effects of the large- 
scale field on the experimental exposure. 

Due to these effects, some regions of the sky have increased 
probability to contribute to the cosmic rays observed on Earth, 
and some other are disfavored. The non-uniform character of 
the mapping strengthens in the case of the heavy primary com- 
position. There exist directions, for which the Earth detectors 
are blind, i.e. even if there is a strong UHECR source in such 
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Fig. 5. The sky distribution of the log 10 of the ratio of the number of events 
at the Galaxy border to the one on Earth, for the HMR bisymmetric even parity 
model [16]. The gray color represents the regions on the sky, to which the 
Earth detectors are (almost) blind, like the one towards the Galactic Center 
on the lower plot. The arrival direction distribution on Earth is isotropic. The 
assumed cosmic ray primary type is indicated above each map. The thick 
dashed line denotes the supergalactic plane. 



a "discriminated" sky region, heavy nuclei from this source 
will not reach the Earth. The highly non-uniform mapping in 
the case of pure iron primary composition is demonstrated by 
the lower plots on Figs. [5] [6] and [7] If the primary cosmic 
rays above 40 EeV are iron nuclei, and the GMF structure is 
adequately described by one of the assumed large-scale GMF 
models, the regions that would effectively contribute to the 
cosmic ray flux on Earth represent only a small fraction of 
the 47r solid angle. 

A. Correlation scan using backtracked directions 

To quantify these large-scale GMF effects on the mapping 
between the arrival directions on Earth and those at the 
Galaxy border in the case of the Pierre Auger Observatory, 
we performed a correlation analysis using backtracked arrival 
directions of simulated events described in Sec. [TV) and the 
694 AGN at redshift z max < 0.024 from the VCV catalogue. 
The simulated events have been divided sequentially into 
samples with the same number of events (81) above 40 EeV as 
in the Auger data. The employed set and ranges of parameters 
were also identical to the ones from [11], [12]. The lower 
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Fig. 6. Idem as Fig. [5] for the HMR axisymmetric odd parity model [16]. Fig. 7. Idem as Fig.[5]for the PS model [19] version by KachelrieB et al. [20]. 



threshold energy giving maximal correlation was scanned 
within a range > 40 EeV. The maximum redshift z max 
scan was performed within a range < z max < 0.024, in 
steps of 0.001. The scan in maximum angular distance d max 
was made within a range 1° < <i max < 8°, in steps of 0.1°. 

We will focus here on the minimum probability values, 
obtained during the scan. The c.d.f. of decimal logarithm 
of the corresponding cumulative binomial probability P m { n of 
reaching this level of correlation under isotropy for the HMR 
BSS_S model are shown on Fig. [8] The level of the minimum 
probability reached in the Auger data is also indicated. The 
results of the correlation scan for the three assumed large- 
scale GMF models and four primary mass compositions are 
summarized in Table HH 

Since for the different GMF model/primary mass assump- 
tions the backtracked directions correlate with the catalogue 
objects in particular privileged regions in the sky, our scan 
results depend strongly on those assumptions. Under some 
configurations, the backtracked directions lie near the super- 
galactic plane, where the catalog object density is higher than 
in average, which increases the correlation probability. For the 
PS model [19] version by KachelrieB et al [20], the scanned 
probability minimum is significantly less deep than for the two 
other spiral-field-only models for any assumed primary mass 
composition, except for protons, where one obtains nearly 
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Fig. 8. Scan results for the simulated event samples, backtracked under the 
HMR bisymmetric even type model [16]. The four cumulative distributions 
of lo 

§io ^min correspond to primary assumptions of pure protons, carbon, 
silicon, or iron nuclei. The vertical dashed line indicates the value obtained 
from the scan of the Auger events [11], [12] on Earth. 



the same level of correlation for all models. The model with 
additional halo field components is clearly less compatible 
with the observed correlation with the nearby AGN [11], [12], 
unless the UHECR flux contribution from these objects (or 
objects with similar spatial distribution) is highly non-uniform 
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HMR bisymmetric even parity model 



p 


1 


1148 


-4.20 


-2.58 


-1.47 


c 


6 


1148 


-6.26 


-3.82 


-2.34 


Si 


14 


1147 


-10.49 


-7.03 


-4.53 


Fe 


26 


1148 


-8.10 


-5.22 


-3.13 



HMR axisymmetric odd parity model 



p 


1 


1148 


-4.48 


-2.70 


-1.57 


c 


6 


1148 


-8.69 


-5.78 


-3.70 


Si 


14 


1148 


-12.56 


-8.50 


-5.30 


Fe 


26 


1148 


-8.98 


-5.89 


-3.85 




PS model version by Kachelriefi et al. 


P 


1 


308 


-4.46 


-2.66 


-1.45 


C 


6 


308 


-4.55 


-2.55 


-1.43 


Si 


14 


308 


-5.21 


-2.92 


-1.72 


Fe 


26 


308 


-4.46 


-2.52 


-1.31 



TABLE II 

Values of log 10 Pmin at the indicated percentiles of the c.d.f., 

FOR THE ASSUMED GMF MODELS AND PRIMARY COMPOSITIONS. 



(for example, if there are only few powerful sources, one of 
which is located by chance in one of those regions in the sky, 
where the large-scale field "enhances" the Auger exposure). 

The presented approach has therefore a potential of pro- 
viding constraints that can be used to discriminate between 
the GMF models and/or the primary charge. It is now being 
complemented by a forward-tracking of cosmic rays to the 
Earth, for a number of plausible UHECR sources scenarios. 

VII. Conclusions 

We have studied the propagation of highest energy 
cosmic rays in the Galactic magnetic field, for different 
assumptions on the large-scale field structure and/or primary 
cosmic ray mass composition. The knowledge of the field 
distribution in the Milky Way is limited by the sensitivity 
and angular resolution of the present-day instruments. The 
Pierre Auger Observatory data provide the measurements 
of the primary cosmic ray properties with unprecedented 
level of collected statistics and reconstruction accuracy. This 
allows a complementary and independent way of probing the 
Galactic field structure. The observed correlation between 
the highest energy Auger events and the nearby AGN may 
give additional hints about the GMF amplitude and orientation. 

We have investigated deflection patterns of cosmic rays 
above 40EeV, underwent during their propagation in the 
Galactic field under three distinctively different GMF models. 
The exact deflection value in the regular component of 
the field depends strongly on the arrival direction on Earth 
of a cosmic ray, with the corresponding position angle of 
deflection differing from one assumed field distribution to 
another. In addition, the spectrometric character of deflections 
of the UHECR in the large-scale Galactic magnetic field 
gives rise to the aligned structures of events coming from 
a UHECR source (thread-like multiplets), that can be used 



to discriminate between GMF models (see [30] and the 
references therein). 



Magnetic deflections and (de-)focusing effects of the Galac- 
tic magnetic field modify the exposure of an experiment 
to the extragalactic sky. Our studies show that the UHECR 
picture observed on Earth is sensitive to the Galactic field 
distribution and the primary cosmic ray composition. Though 
the reconstruction of the field is easier in the case of light 
primary mass composition, in the case of heavy nuclei the 
lensing effects of the Galactic field on the exposure bring 
stronger constraints on the list of potential UHECR source 
candidates. The presented analysis of the correspondence 
between the arrival direction distributions on Earth and at 
the Galactic border, and of matching of the latter with the 
VCV catalogue objects will be complemented by the forward- 
tracking of cosmic rays for a number of plausible UHECR 
sources scenarios. This will allow for direct comparison with 
the observed AGN correlation. 
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